function out = yieldCurveDecom_threshold(setup,output,theta2,numSim)
% Information: we use the value of sigma for which theta11 is estimated to compute the bond yields.
% To simulate the factors, we use the value of sigma from the SEVAR model
% The conditional expectations of the short rate are computed by simulation
tmp                = output.model;
nm                 = setup.nm;
setupAll           = setup;
setupAll.matSelect = 1:max(setup.matSelect);
if output.model.modelNum == 1
    [modelAll,errorMes]= solveQTSM(tmp.alpha,tmp.beta,tmp.psi,tmp.mu,tmp.phi,tmp.sigmaxtil,setupAll);
elseif output.model.modelNum == 2
    [modelAll,errorMes] = solveQTSM_constRsmooth(tmp.alpha,tmp.rhor,tmp.beta,tmp.psi,tmp.mu,tmp.phi,tmp.sigmaxtil,setup);
end
modelAll.modelNum = output.model.modelNum;
modelAll.factorSignResOn = 0;
modelAll.macros   = output.macros;
modelAll.rLag     = output.model.rLag;
modelAll.rhor     = output.model.rhor;
econAct           = setup.econAct;

% Getting the P dynamics. Here we use the estimates from the VAR model in step 3
factors          = [output.macros'; output.factors];
[nx,T]           = size(factors);
[h0_1,h0_2,hx_1,hx_2,sigmaVar,gama0,gamaz,gamax,sigmaz] = unfoldTheta2_threshold_macro(theta2,nx);
maxMat           = max(modelAll.matSelect);
%% The distribution of the factors under P 
seedNum = 1;
rng(seedNum,'twister')
factorSim = nan(nx,numSim,maxMat);
tmp       = nan(numSim,maxMat);
zSim      = nan(numSim,maxMat);
rExp      = nan(maxMat,T);
epsilonx                     = randn(nx,numSim,maxMat);
epsilonx(:,numSim/2+1:end,:) = -epsilonx(:,1:numSim/2,:);
epsilonz                     = randn(numSim,maxMat);
epsilonz(numSim/2+1:end,:)   = -epsilonz(1:numSim/2,:);
disp('Computing conditional expected short rates by Monte Carlo integration .... ')
vectorOn                     = 1;
for t = 1:T
    if mod(t,50) == 0
        disp(['t = ',num2str(t)])
    end
    % Vectorized version
    if vectorOn == 1
        for i = 1:maxMat
            if i == 1
                factorSim(:,:,i) = repmat(factors(:,t),1,numSim);
                zSim(:,i)        = econAct(t,1);
            else
                boom  = (zSim(:,i-1) >= 0)';
                reces = (zSim(:,i-1) < 0)';
                factorSim(:,:,i) = ...
                    repmat(boom,nx,1).*repmat(h0_1,1,numSim)...
                    + repmat(reces,nx,1).*repmat(h0_2,1,numSim) ...
                    + repmat(boom,nx,1).*(hx_1*factorSim(:,:,i-1)) ...
                    + repmat(reces,nx,1).*(hx_2*factorSim(:,:,i-1))...
                    + sigmaVar*epsilonx(:,:,i); %randn(nx,1);
                zSim(:,i) = repmat(gama0,numSim,1) + gamaz*zSim(:,i-1) ...
                    + (gamax'*factorSim(:,:,i-1))' + sigmaz*epsilonz(:,i); %randn(1,1);
            end
            if nx == 3
                xxTmp = [factorSim(1,:,i).*factorSim(:,:,i); ...
                         factorSim(2,:,i).*factorSim(:,:,i); ...
                         factorSim(3,:,i).*factorSim(:,:,i)];
            elseif nx == 4
                xxTmp = [factorSim(1,:,i).*factorSim(:,:,i); ...
                         factorSim(2,:,i).*factorSim(:,:,i); ...
                         factorSim(3,:,i).*factorSim(:,:,i); ...
                         factorSim(4,:,i).*factorSim(:,:,i)];
            elseif nx == 5
                xxTmp = [factorSim(1,:,i).*factorSim(:,:,i); ...
                         factorSim(2,:,i).*factorSim(:,:,i); ...
                         factorSim(3,:,i).*factorSim(:,:,i); ...
                         factorSim(4,:,i).*factorSim(:,:,i); ...
                         factorSim(5,:,i).*factorSim(:,:,i)];
            else
                error('in yieldCurveDecom_threshold, nx > 5');
            end
            tmp(:,i) = modelAll.r0 + modelAll.rx*factorSim(:,:,i) + modelAll.rxx*xxTmp;
        end
    else
        for s=1:numSim
            for i = 1:maxMat
                if i == 1
                    factorSim(:,s,1) = factors(:,t);
                    zSim(s,1)        = econAct(t,1);
                else
                    if zSim(s,i-1) >= 0
                        factorSim(:,s,i) = h0_1 + hx_1*factorSim(:,s,i-1) + sigmaVar*epsilonx(:,s,i); %randn(nx,1);
                    else
                        factorSim(:,s,i) = h0_2 + hx_2*factorSim(:,s,i-1) + sigmaVar*epsilonx(:,s,i); %randn(nx,1);
                    end
                    zSim(s,i) = gama0 + gamaz*zSim(s,i-1) + gamax'*factorSim(:,s,i-1) + sigmaz*epsilonz(s,i); %randn(1,1);
                end
                tmp(s,i) = modelAll.r0 + modelAll.rx*factorSim(:,s,i) ...
                    + modelAll.rxx*kron3(factorSim(:,s,i),factorSim(:,s,i));
            end
        end
    end
    
    if setup.modelNum == 1
        rExp(:,t) = mean(tmp,1)';
    elseif setup.modelNum == 2
        rExp(1,t) = modelAll.rhor*modelAll.rLag(t,1)+mean(tmp(:,1),1);
        for i=2:maxMat
            rExp(i,t) = modelAll.rhor*rExp(i-1,t)+ mean(tmp(:,i),1);
        end
    end
end

%% Fitted yields and term premia
% tp(n)_t = y(n)_t - 1/n*sum(E_t[r_t+i],i=0,1,...n-1)
ny         = size(modelAll.matSelect,2);
yHat       = nan(ny,T);
termPrem   = nan(ny,T);
rExpAvg    = nan(ny,T);
for t=1:T
    yHat(:,t) = modelFit(modelAll,factors(nm+1:end,t),t);
    for i=1:ny
        rExpAvg(i,t) = mean(rExp(1:modelAll.matSelect(1,i),t));
        termPrem(i,t) = yHat(i,t) - rExpAvg(i,t);
    end
end

%% Convert into forwards
fHat    = nan(ny,T);
fwdPrem = nan(ny,T);
for i = 1:ny-1
    fHat(i,:)    = (i+1)*yHat(i+1,:) - i*yHat(i,:);
    fwdPrem(i,:) = (i+1)*termPrem(i+1,:) - i*termPrem(i,:);
end

%% Output
out.yieldsHat     = yHat;
out.termPremia    = termPrem*10000;
out.rExp          = rExp;
out.forwardsHat   = fHat;
out.forwardPermia = fwdPrem*10000;
out.matSelect     = modelAll.matSelect;
out.numSim        = numSim;

end